HS project

Figures

This file is used to prepare the figures for the paper.

library(dplyr)
library(patchwork)
library(ggplot2)
library(ComplexHeatmap)
library(org.Mm.eg.db)

.libPaths()

## [1] "/usr/local/lib/R/library"

Preparation

Here are the folders where analyzes are stored :

data_dir = "./.."
list.files(data_dir)

## [1] "0_intro"      "1_metadata"   "2_individual" "3_combined"   "4_zoom"      
## [6] "5_figures"    "LICENSE"      "index.html"   "index_layout"

We load the dataset containing all cells :

sobj = readRDS(paste0(data_dir, "/3_combined/hs_hd_sobj.rds"))
sobj

## An object of class Seurat 
## 20003 features across 12111 samples within 1 assay 
## Active assay: RNA (20003 features, 2000 variable features)
##  6 dimensional reductions calculated: RNA_pca, RNA_pca_38_tsne, RNA_pca_38_umap, harmony, harmony_38_umap, harmony_38_tsne

These are all the samples analyzed :

sample_info = readRDS(paste0(data_dir, "/1_metadata/hs_hd_sample_info.rds"))

# Nb cells by dataset
to_plot = table(sobj$sample_identifier) %>%
  as.data.frame.table(., stringsAsFactors = FALSE) %>%
  `colnames<-`(c("sample_identifier", "nb_cells")) %>%
  dplyr::left_join(x = ., y = sample_info, by = "sample_identifier") 

# patchwork
plot_list = aquarius::fig_plot_gb(to_plot, title = "Available datasets")
patchwork::wrap_plots(plot_list) +
  patchwork::plot_layout(design = "A\nB", heights = c(0.1,5)) &
  ggplot2::theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 15))

These are the custom colors for cell populations :

color_markers = readRDS(paste0(data_dir, "/1_metadata/hs_hd_color_markers.rds"))
color_markers = color_markers[names(color_markers) != "melanocytes"]
ors_color = color_markers["ORS"]
color_markers["ORS"] = color_markers["IFE"] 
color_markers["IFE"] = ors_color
color_markers["B cells"] = "chocolate3"
rm(ors_color)

data.frame(cell_type = names(color_markers),
           color = unlist(color_markers)) %>%
  ggplot2::ggplot(., aes(x = cell_type, y = 0, fill = cell_type)) +
  ggplot2::geom_point(pch = 21, size = 5) +
  ggplot2::scale_fill_manual(values = unlist(color_markers), breaks = names(color_markers)) +
  ggplot2::theme_classic() +
  ggplot2::theme(legend.position = "none",
                 axis.line = element_blank(),
                 axis.title = element_blank(),
                 axis.ticks = element_blank(),
                 axis.text.y = element_blank(),
                 axis.text.x = element_text(hjust = 1, angle = 20))

We define custom colors for sample type :

sample_type_colors = setNames(nm = levels(sample_info$sample_type),
                              c("#C55F40", "#2C78E6"))

data.frame(cell_type = names(sample_type_colors),
           color = unlist(sample_type_colors)) %>%
  ggplot2::ggplot(., aes(x = cell_type, y = 0, fill = cell_type)) +
  ggplot2::geom_point(pch = 21, size = 5) +
  ggplot2::scale_fill_manual(values = unlist(sample_type_colors), breaks = names(sample_type_colors)) +
  ggplot2::theme_classic() +
  ggplot2::theme(legend.position = "none",
                 axis.line = element_blank(),
                 axis.title = element_blank(),
                 axis.ticks = element_blank(),
                 axis.text.y = element_blank())

We set a background color :

bg_color = "gray94"

This is the correspondence between cell types and cell families, and custom colors to color cells by cell family :

custom_order_cell_type = data.frame(
  cell_type = names(color_markers),
  cell_family = c(rep("immune cells", 5),
                  rep("matrix", 5),
                  rep("non matrix", 5)),
  stringsAsFactors = FALSE)
custom_order_cell_type$cell_type = factor(custom_order_cell_type$cell_type,
                                          levels = custom_order_cell_type$cell_type)
rownames(custom_order_cell_type) = custom_order_cell_type$cell_type

family_color = c("immune cells" = "slateblue1",
                 "matrix" = "mediumseagreen",
                 "non matrix" = "firebrick3")

We load markers to display on a dotplot to assess cell type annotation :

dotplot_markers = readRDS(paste0(data_dir, "/1_metadata/hs_hd_dotplot_markers.rds"))
dotplot_markers = dotplot_markers[names(dotplot_markers) != "melanocytes"]
lengths(dotplot_markers)

##      CD4 T cells      CD8 T cells Langerhans cells      macrophages 
##                2                2                2                2 
##          B cells          cuticle           cortex          medulla 
##                2                2                2                2 
##              IRS    proliferative              IBL              ORS 
##                2                2                2                2 
##              IFE             HFSC        sebocytes 
##                2                2                2

Custom functions to display gene expression on the heatmap :

color_fun = function(one_gene) {
  gene_range = range(ht_annot[, one_gene])
  gene_palette = circlize::colorRamp2(colors = c("#FFFFFF", aquarius::color_gene[-1]),
                                      breaks = seq(from = gene_range[1], to = gene_range[2],
                                                   length.out = length(aquarius::color_gene)))
  return(gene_palette)
}

All samples

Settings

This is the projection name to visualize cells :

name2D = "harmony_38_tsne"
name2D_atlas = name2D

Preparation

We make a low resolutive clustering for the heatmap :

sobj = Seurat::FindClusters(sobj, resolution = 0.4)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 12111
## Number of edges: 475472
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9421
## Number of communities: 17
## Elapsed time: 1 seconds

length(levels(sobj$seurat_clusters))

## [1] 17

We define cluster type and cluster family :

cluster_type = table(sobj$cell_type, sobj$seurat_clusters) %>%
  prop.table(., margin = 2) %>%
  apply(., 2, which.max)
cluster_type = setNames(nm = names(cluster_type),
                        levels(sobj$cell_type)[cluster_type])

sobj$cluster_type = cluster_type[sobj$seurat_clusters]
sobj$cluster_type = factor(sobj$cluster_type,
                           levels = levels(sobj$cell_type))
sobj$cluster_family = custom_order_cell_type[sobj$cluster_type, "cell_family"]
sobj$cluster_family = factor(sobj$cluster_family,
                             levels = names(family_color))

Figures

Project name :

# Random order
set.seed(1234)
rnd_order = sample(colnames(sobj), replace = FALSE, size = ncol(sobj))

# Extract coordinates
cells_coord = sobj@reductions[[name2D]]@cell.embeddings %>%
  as.data.frame() %>%
  `colnames<-`(c("Dim1", "Dim2"))
cells_coord$project_name = sobj$project_name
cells_coord = cells_coord[(rnd_order), ]

# Plot
ggplot2::ggplot(cells_coord, aes(x = Dim1, y = Dim2, col = project_name)) +
  ggplot2::geom_point(size = 0.5) +
  ggplot2::scale_color_manual(values = sample_info$color,
                              breaks = sample_info$project_name) +
  ggplot2::theme_void() +
  ggplot2::theme(aspect.ratio = 1,
                 legend.position = "none")

Cluster :

Seurat::DimPlot(sobj, reduction = name2D, pt.size = 0.4,
                group.by = "seurat_clusters", label = TRUE) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() +
  Seurat::NoLegend()

Cell type annotation :

Seurat::DimPlot(sobj, reduction = name2D, pt.size = 0.5,
                group.by = "cell_type", cols = color_markers) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() +
  Seurat::NoLegend()

Cluster family annotation :

Seurat::DimPlot(sobj, reduction = name2D, pt.size = 0.5,
                group.by = "cluster_family", cols = family_color) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() +
  Seurat::NoLegend()

Cell type annotation split by condition :

plot_list = aquarius::plot_split_dimred(sobj, reduction = name2D,
                                        group_by = "cell_type",
                                        group_color = color_markers,
                                        split_by = "sample_type",
                                        split_color = sample_type_colors,
                                        bg_color = bg_color)

patchwork::wrap_plots(plot_list) &
  Seurat::NoLegend()

Gene expression to assess annotation :

genes = c("PTPRC", "MSX2", "KRT14")
names(genes) = c("immune cells", "matrix cells", "non-matrix cells")

plot_list = lapply(c(1:length(genes)), FUN = function(gene_id) {
  gene = genes[[gene_id]]
  pop = names(genes)[gene_id]
  
  sobj$my_gene = Seurat::FetchData(sobj, gene)[, 1] %>%
    aquarius::run_rescale(., new_min = 0, new_max = 10)
  
  Seurat::FeaturePlot(sobj, features = "my_gene", reduction = name2D) +
    ggplot2::scale_color_gradientn(colors = aquarius::color_gene,
                                   breaks = seq(0, 10, by = 2.5),
                                   labels = c("min", rep("", 3), "max")) +
    ggplot2::labs(title = gene) + 
    # subtitle = pop) +
    ggplot2::theme(aspect.ratio = 1,
                   plot.title = element_text(hjust = 0.5, size = 17),
                   plot.subtitle = element_text(hjust = 0.5, size = 15),
                   legend.text = element_text(size = 15),
                   legend.position = "none") +
    Seurat::NoAxes()
})

plot_list

## [[1]]

## 
## [[2]]

## 
## [[3]]

Barplot by cluster family :

quantif = table(sobj$sample_identifier) %>%
  as.data.frame.table() %>%
  `colnames<-`(c("Sample", "nb_cells"))

aquarius::plot_barplot(df = table(sobj$sample_identifier,
                                  sobj$cluster_family) %>%
                         as.data.frame.table() %>%
                         `colnames<-`(c("Sample", "Cell Type", "Number")),
                       x = "Sample", y = "Number", fill = "Cell Type",
                       position = position_fill()) +
  ggplot2::geom_label(data = quantif, inherit.aes = FALSE,
                      aes(x = .data$Sample, y = 1.05, label = .data$nb_cells),
                      label.size = 0, size = 5) +
  ggplot2::scale_fill_manual(values = unlist(family_color),
                             breaks = names(family_color),
                             name = "Cell Family") +
  ggplot2::scale_y_continuous(breaks = seq(0, 1, by = 0.25),
                              labels = paste0(seq(0, 100, by = 25), sep = " %"),
                              expand = ggplot2::expansion(add = c(0, 0.05))) +
  ggplot2::theme(axis.title.y = element_blank(),
                 axis.line.x = element_line(colour = "lightgray"),
                 text = element_text(size = 15),
                 legend.position = "none")

Heatmap of cluster proportion by sample :

group_by = "seurat_clusters"

cluster_by_sample = table(sobj$sample_identifier,
                          sobj@meta.data[, group_by]) %>%
  prop.table(margin = 1) %>%
  as.matrix()

## Representative markers
cluster_markers = c("PTPRC", "CD3E", "AIF1", "MSX2", "DIO2", "GPX2", "KRT16")
ht_annot = Seurat::GetAssayData(sobj, assay = "RNA", slot = "data")[cluster_markers, ] %>%
  Matrix::t() %>% as.data.frame()
ht_annot$clusters = sobj@meta.data[, group_by]
ht_annot = ht_annot %>%
  dplyr::group_by(clusters) %>%
  dplyr::summarise_all(funs('mean' = mean)) %>%
  as.data.frame() %>%
  dplyr::select(-clusters) %>%
  `colnames<-`(c(cluster_markers))

## Bottom annotation : average gene expression
ha_bottom = ComplexHeatmap::HeatmapAnnotation(df = ht_annot,
                                              which = "column",
                                              show_legend = FALSE,
                                              col = setNames(nm = cluster_markers,
                                                             lapply(cluster_markers, FUN = color_fun)),
                                              annotation_name_side = "right")

## Right annotation : number of cells by dataset
ht_annot = table(sobj$sample_identifier) %>%
  as.data.frame.table() %>%
  `colnames<-`(c("sample_identifier", "nb_cells")) %>%
  `rownames<-`(.$sample_identifier) %>%
  dplyr::select(-sample_identifier)

ha_right = ComplexHeatmap::HeatmapAnnotation(
  df = ht_annot,
  which = "row",
  show_legend = TRUE,
  annotation_name_side = "top",
  col = list(nb_cells  = circlize::colorRamp2(colors = RColorBrewer::brewer.pal(name = "Greys", n = 9),
                                              breaks = seq(from = range(ht_annot$nb_cells)[1],
                                                           to = range(ht_annot$nb_cells)[2],
                                                           length.out = 9))))

## Left annotation : gender
ha_left = ComplexHeatmap::HeatmapAnnotation(
  gender = sample_info$gender,
  which = "row",
  show_legend = TRUE,
  annotation_name_side = "bottom",
  col = list(gender = setNames(nm = c("F", "M"),
                               c("lightcyan3", "navyblue"))))

## Top annotation : main cell type in this cluster
ht_annot = table(sobj$sample_identifier,
                 sobj@meta.data[, group_by]) %>%
  prop.table(margin = 1) %>%
  as.matrix()

ht_annot = table(sobj$cell_type,
                 sobj@meta.data[, group_by]) %>%
  prop.table(., margin = 2) %>%
  apply(., 2, which.max)
ht_annot = data.frame(row.names = names(ht_annot),
                      cell_type = names(color_markers)[ht_annot],
                      stringsAsFactors = FALSE)
ht_annot = dplyr::left_join(ht_annot, custom_order_cell_type, by = "cell_type") %>%
  # Simplification for matrix
  dplyr::mutate(cell_type = ifelse(cell_type %in% c("medulla", "cortex", "cuticle"), yes = "hair shaft", no = cell_type)) %>%
  # Simplification for T cells
  dplyr::mutate(cell_type = ifelse(cell_type %in% c("CD4 T cells", "CD8 T cells"), yes = "T cells", no = cell_type)) %>%
  # Simplification for APC
  dplyr::mutate(cell_type = ifelse(cell_type %in% c("Langerhans cells", "macrophages"), yes = "APC", no = cell_type)) %>%
  # Add color
  dplyr::mutate(color = as.character(color_markers[cell_type])) %>%
  dplyr::mutate(color = ifelse(cell_type == "hair shaft", yes = "#FFB6C1", no = color)) %>%
  dplyr::mutate(color = ifelse(cell_type == "T cells", yes = "#8A6EE6", no = color)) %>%
  dplyr::mutate(color = ifelse(cell_type == "APC", yes = "#9CAA4B", no = color))

ha_top = ComplexHeatmap::HeatmapAnnotation(
  cell_type = ht_annot$cell_type,
  # cell_family = ht_annot$cell_family,
  which = "column",
  show_legend = TRUE,
  annotation_name_side = "left",
  col = list(cell_type = setNames(nm = ht_annot$cell_type,
                                  ht_annot$color)
             #, cell_family = family_color
  ))

## Assemble heatmap
ht = ComplexHeatmap::Heatmap(cluster_by_sample,
                             heatmap_legend_param = list(title = "Proportion",
                                                         col = c("#2166AC", "#F7F7F7", "#B2182B")),
                             # bottom_annotation = ha_bottom,
                             right_annotation = ha_right,
                             left_annotation = ha_left,
                             top_annotation = ha_top,
                             cluster_rows = TRUE,
                             cluster_columns = TRUE,
                             row_title = "Sample",
                             row_names_gp = grid::gpar(names = sample_info$sample_identifier,
                                                       col = sample_info$color,
                                                       fontface = "bold"),
                             column_title = "Cluster",
                             column_names_centered = TRUE,
                             row_names_side = "left",
                             column_names_side = "top",
                             column_names_rot = 0)

## Draw !
ComplexHeatmap::draw(ht, merge_legends = TRUE)

Dotplot :

plot_list = aquarius::plot_dotplot(sobj,
                                   markers = c("PTPRC",
                                               "CD3E", "CD4",
                                               "CD3E", "CD8A",
                                               "CD207", "CPVL", 
                                               "TREM2", "MSR1",
                                               "CD79A", "CD79B",
                                               "KRT32", "KRT35",
                                               "KRT31", "PRR9",
                                               "BAMBI", "ADLH1A3",
                                               "KRT71", "KRT73",
                                               "TOP2A", "MCM5",
                                               "KRT16", "KRT6C",
                                               "KRT15", "GPX2",
                                               "SPINK5", "LY6D",
                                               "DIO2", "TCEAL2",
                                               "CLMP", "PPARG"),
                                   assay = "RNA", column_name = "cell_type", nb_hline = 0) +
  ggplot2::scale_color_gradientn(colors = aquarius::color_gene) +
  ggplot2::theme(legend.position = "right",
                 legend.box = "vertical",
                 axis.title = element_blank(),
                 axis.ticks.x = element_blank(),
                 axis.text.x = element_blank(),
                 axis.line.x = element_blank(),
                 plot.margin = unit(rep(0, 4), "cm"))

p = ggplot2::ggplot(custom_order_cell_type, aes(x = cell_type, y = 0)) +
  ggplot2::geom_point(size = 0) +
  ggplot2::geom_segment(aes(x = 0.5, xend = 5.5, y = 0, yend = 0), size = 6, col = family_color["immune cells"]) +
  ggplot2::geom_segment(aes(x = 5.5, xend = 10.5, y = 0, yend = 0), size = 6, col = family_color["matrix"]) +
  ggplot2::geom_segment(aes(x = 10.5, xend = 15.5, y = 0, yend = 0), size = 6, col = family_color["non matrix"]) +
  ggplot2::scale_y_continuous(expand = c(0,0), limits = c(0,0)) +
  ggplot2::theme_classic() +
  ggplot2::theme(axis.text.y = element_blank(),
                 axis.ticks.y = element_blank(),
                 axis.title = element_blank(),
                 axis.line.y = element_blank(),
                 axis.text.x = element_text(angle = 45, hjust = 1, size = 10, color = "black"),
                 plot.margin = unit(rep(0, 4), "cm"))

plot_list = patchwork::wrap_plots(plot_list, p,
                                  ncol = 1, heights = c(25, 1))
plot_list

Immune cells

Settings

We load the immune cells dataset :

sobj_ic = readRDS(paste0(data_dir, "/4_zoom/1_zoom_immune/immune_cells_sobj.rds"))
sobj_ic

## An object of class Seurat 
## 15121 features across 2329 samples within 1 assay 
## Active assay: RNA (15121 features, 2000 variable features)
##  6 dimensional reductions calculated: RNA_pca, RNA_pca_20_tsne, RNA_pca_20_umap, harmony, harmony_20_umap, harmony_20_tsne

This is the projection name to visualize cells :

name2D = "harmony_20_tsne"

To represent results from differential expression, we load the analyses results :

list_results = readRDS(paste0(data_dir, "/4_zoom/1_zoom_immune/immune_cells_list_results.rds"))

lapply(list_results, FUN = names)

## $`Langerhans cells`
## [1] "mark" "gsea"
## 
## $macrophages
## [1] "mark"       "enrichr_hs" "enrichr_hd" "gsea"      
## 
## $`CD4 T cells`
## [1] "mark"       "enrichr_hs" "enrichr_hd" "gsea"      
## 
## $`CD8 T cells`
## [1] "mark"       "enrichr_hs" "enrichr_hd" "gsea"

Preparation

We defined cluster type and cluster family :

cluster_type = table(sobj_ic$cell_type, sobj_ic$seurat_clusters) %>%
  prop.table(., margin = 2) %>%
  apply(., 2, which.max)
cluster_type = setNames(nm = names(cluster_type),
                        levels(sobj_ic$cell_type)[cluster_type])

sobj_ic$cluster_type = cluster_type[sobj_ic$seurat_clusters]
sobj_ic$cluster_type = factor(sobj_ic$cluster_type,
                              levels = levels(sobj_ic$cell_type))
sobj_ic$cluster_family = custom_order_cell_type[sobj_ic$cluster_type, "cell_family"]
sobj_ic$cluster_family = factor(sobj_ic$cluster_family,
                                levels = names(family_color))

Figures

Control cells on the full atlas :

sobj$is_immune = (colnames(sobj) %in% colnames(sobj_ic))

Seurat::DimPlot(sobj, reduction = name2D_atlas, pt.size = 0.000001,
                group.by = "is_immune", order = "TRUE") +
  ggplot2::scale_color_manual(values = c(family_color[["immune cells"]], bg_color),
                              breaks = c(TRUE, FALSE)) +
  ggplot2::labs(title = "Immune cells",
                subtitle = paste0(ncol(sobj_ic), " cells")) +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_text(hjust = 0.5, face = "bold"),
                 plot.subtitle = element_text(hjust = 0.5)) +
  Seurat::NoAxes() + Seurat::NoLegend()

Violin plot of IL1B in macrophages :

subsobj = subset(sobj_ic, seurat_clusters == 2)
table(subsobj$sample_type)

## 
##  HS  HD 
## 378  31

il1b_hs = subsobj@assays$RNA@data["IL1B", subsobj$sample_type == "HS"]
il1b_hd = subsobj@assays$RNA@data["IL1B", subsobj$sample_type == "HD"]
il1b_hs_VS_il1b_hd = stats::t.test(il1b_hs, il1b_hd)
il1b_hs_VS_il1b_hd

## 
##  Welch Two Sample t-test
## 
## data:  il1b_hs and il1b_hd
## t = 2.3206, df = 35.801, p-value = 0.02612
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.05066074 0.75429252
## sample estimates:
## mean of x mean of y 
## 0.9256690 0.5231924

Seurat::VlnPlot(subsobj, group.by = "sample_type",
                features = "IL1B", cols = sample_type_colors) +
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

Split by sample :

Seurat::VlnPlot(subsobj, group.by = "sample_identifier",
                features = "IL1B", cols = sample_info$color) +
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

Violin plot of IL1B in macrophages :

il6_hs = subsobj@assays$RNA@data["IL6", subsobj$sample_type == "HS"]
il6_hd = subsobj@assays$RNA@data["IL6", subsobj$sample_type == "HD"]
il6_hs_VS_il6_hd = stats::t.test(il6_hs, il6_hd)
il6_hs_VS_il6_hd

## 
##  Welch Two Sample t-test
## 
## data:  il6_hs and il6_hd
## t = 2.3591, df = 377, p-value = 0.01883
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.001453521 0.016004661
## sample estimates:
##   mean of x   mean of y 
## 0.008729091 0.000000000

Seurat::VlnPlot(subsobj, group.by = "sample_type",
                features = "IL6", cols = sample_type_colors) +
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

Violin plot of IL1B in macrophages :

subsobj = subset(sobj_ic, seurat_clusters == 2)
table(subsobj$sample_type)

## 
##  HS  HD 
## 378  31

tnf_hs = subsobj@assays$RNA@data["TNF", subsobj$sample_type == "HS"]
tnf_hd = subsobj@assays$RNA@data["TNF", subsobj$sample_type == "HD"]
tnf_hs_VS_tnf_hd = stats::t.test(tnf_hs, tnf_hd)
tnf_hs_VS_tnf_hd

## 
##  Welch Two Sample t-test
## 
## data:  tnf_hs and tnf_hd
## t = 3.8395, df = 45.546, p-value = 0.0003786
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.1140962 0.3657054
## sample estimates:
##  mean of x  mean of y 
## 0.31784960 0.07794879

Seurat::VlnPlot(subsobj, group.by = "sample_type",
                features = "TNF", cols = sample_type_colors) +
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

Split by sample :

Seurat::VlnPlot(subsobj, group.by = "sample_identifier",
                features = "TNF", cols = sample_info$color) +
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

Violin plot of GZMA in CD4 T cells :

subsobj = subset(sobj_ic, seurat_clusters %in% c(0,10))
table(subsobj$sample_type)

## 
##  HS  HD 
## 569  90

gzma_hs = subsobj@assays$RNA@data["GZMA", subsobj$sample_type == "HS"]
gzma_hd = subsobj@assays$RNA@data["GZMA", subsobj$sample_type == "HD"]
gzma_hs_VS_gzma_hd = stats::t.test(gzma_hs, gzma_hd)
gzma_hs_VS_gzma_hd

## 
##  Welch Two Sample t-test
## 
## data:  gzma_hs and gzma_hd
## t = 14.755, df = 172.51, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.225578 1.604097
## sample estimates:
## mean of x mean of y 
## 1.8456108 0.4307735

Seurat::VlnPlot(subsobj, group.by = "sample_type",
                features = "GZMA", cols = sample_type_colors) +
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

We represent some genes split by sample type :

plot_list = lapply(c("IL1B", "GZMA", "IFNG", "IL17A"), FUN = function(one_gene) {
  p = aquarius::plot_split_dimred(sobj_ic,
                                  reduction = name2D,
                                  split_by = "sample_type",
                                  color_by = one_gene,
                                  color_palette = c("gray70", "#FDBB84", "#EF6548", "#7F0000", "black"),
                                  main_pt_size = 0.6,
                                  bg_pt_size = 0.6,
                                  order = TRUE,
                                  bg_color = "gray95")
  p = patchwork::wrap_plots(p, nrow = 1) +
    patchwork::plot_layout(guides = "collect") +
    ggplot2::theme(legend.position = "right") &
    ggplot2::theme(plot.subtitle = element_blank())
  return(p)
})

patchwork::wrap_plots(plot_list, ncol = 2)

Barplot by cluster type :

quantif = table(sobj_ic$sample_identifier) %>%
  as.data.frame.table() %>%
  `colnames<-`(c("Sample", "nb_cells"))

aquarius::plot_barplot(df = table(sobj_ic$sample_identifier,
                                  sobj_ic$cluster_type) %>%
                         as.data.frame.table() %>%
                         `colnames<-`(c("Sample", "Cell Type", "Number")),
                       x = "Sample", y = "Number", fill = "Cell Type",
                       position = position_stack()) +
  ggplot2::geom_label(data = quantif, inherit.aes = FALSE,
                      aes(x = .data$Sample, y = 50 + .data$nb_cells, label = .data$nb_cells),
                      label.size = 0, size = 5) +
  ggplot2::scale_fill_manual(values = unlist(color_markers),
                             breaks = names(color_markers),
                             name = "Cell Type") +
  ggplot2::theme(axis.title.y = element_blank(),
                 axis.line.x = element_line(colour = "lightgray"),
                 text = element_text(size = 15),
                 legend.position = "none")

Dotplot of DE genes for macrophages :

subsobj = subset(sobj_ic, cluster_type == "macrophages")
features_oi = c("HLA-DQA2", "HLA-DPA1", "HLA-DRB5",
                "HLA-A", "HLA-C", "B2M",
                "C1QA", "C1QB", "C1QC")

Seurat::DotPlot(subsobj,
                assay = "RNA",
                features = features_oi,
                group.by = "sample_type",
                scale = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::scale_color_gradientn(colors = aquarius::color_gene) +
  ggplot2::theme(axis.title.x = element_blank(),
                 axis.title.y = element_blank(),
                 axis.text.x = element_text(angle = 30, hjust = 1))

Heatmap for macrophages :

# Matrix
mat_expr = Seurat::GetAssayData(subsobj)
mat_expr = mat_expr[features_oi, ]
mat_expr = Matrix::t(mat_expr)
mat_expr = dynutils::scale_quantile(mat_expr) # between 0 and 1
mat_expr = Matrix::t(mat_expr)
dim(mat_expr) # genes x cells

## [1]   9 463

## Colors
list_colors = list()

# Heatmap
list_colors[["expression"]] = rev(RColorBrewer::brewer.pal(name = "RdBu", n = 9))

# Sample annotation (top annotation)
list_colors[["sample_type"]] = sample_type_colors
list_colors[["sample_identifier"]] = setNames(nm = sample_info$sample_identifier,
                                              sample_info$color)
# Cells order
column_order = subsobj@meta.data %>%
  dplyr::arrange(sample_type, sample_identifier) %>%
  rownames()
column_order = match(column_order, rownames(subsobj@meta.data))

# Heatmap
ha_top = HeatmapAnnotation(sample_type = subsobj$sample_type,
                           sample_identifier = subsobj$sample_identifier,
                           col = list(sample_type = list_colors[["sample_type"]],
                                      sample_identifier = list_colors[["sample_identifier"]]))

# Heatmap
ht = Heatmap(as.matrix(mat_expr),
             heatmap_legend_param = list(title = "Expression", at = c(0, 1), 
                                         labels = c("low", "high")),
             col = list_colors[["expression"]],
             top_annotation = ha_top,
             show_column_names = FALSE,
             column_order = column_order,
             column_gap = unit(2, "mm"),
             cluster_rows = FALSE,
             row_title = NULL,
             row_names_gp = grid::gpar(fontsize = 14, fontface = "plain"),
             use_raster = FALSE,
             show_heatmap_legend = TRUE,
             border = TRUE)

ComplexHeatmap::draw(ht,
                     merge_legend = TRUE,
                     heatmap_legend_side = "bottom",
                     annotation_legend_side = "bottom")

Dotplot of DE genes for CD4 T cells :

subsobj = subset(sobj_ic, cluster_type == "CD4 T cells")
# features_oi = c(
#   # HALLMARK_TNFA_SIGNALING_VIA_NFKB
#   "RIPK2", "KLF9", "PER1", "FOS", "BTG1", "PTGER4", "BIRC3", "AREG", "KLF6", "ZFP36", "NFKBIA",
#   # REACTOME_RHO_GTPASE_EFFECTORS
#   "CYBA", "ARPC3", "CFL1", "RPS27", "CALM1", "ACTB",
#   # GOBP_GRANZYME_MEDIATED_PROGRAMMED_CELL_DEATH_SIGNALING_PATHWAY
#   "BNIP3", "GZMA", "GZMB", "LAMP1", "NKG7", "SRGN", "UBE4B",
#   # GOBP_POSITIVE_REGULATION_OF_MEMORY_T_CELL_DIFFERENTIATION
#   "CD46", "HLA-DRA", "HLA-DRB1", "IL12RB1", "IL23A", "IL23R", "TNFSF4",
#   "KLRB1")
features_oi = c("GZMA", "KLRB1", "BTG1", "ZFP36", "NFKBIA", "TXNIP", "CXCR4")

Seurat::DotPlot(subsobj,
                assay = "RNA",
                features = features_oi,
                group.by = "sample_type",
                scale = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::scale_color_gradientn(colors = aquarius::color_gene) +
  ggplot2::theme(axis.title.x = element_blank(),
                 axis.title.y = element_blank(),
                 axis.text.x = element_text(angle = 30, hjust = 1))

Heatmap for CD4 T cells :

# Matrix
mat_expr = Seurat::GetAssayData(subsobj)
mat_expr = mat_expr[features_oi, ]
mat_expr = Matrix::t(mat_expr)
mat_expr = dynutils::scale_quantile(mat_expr) # between 0 and 1
mat_expr = Matrix::t(mat_expr)
dim(mat_expr) # genes x cells

## [1]   7 848

## Colors
list_colors = list()

# Heatmap
list_colors[["expression"]] = rev(RColorBrewer::brewer.pal(name = "RdBu", n = 9))

# Sample annotation (top annotation)
list_colors[["sample_type"]] = sample_type_colors
list_colors[["sample_identifier"]] = setNames(nm = sample_info$sample_identifier,
                                              sample_info$color)
# Cells order
column_order = subsobj@meta.data %>%
  dplyr::arrange(sample_type, sample_identifier) %>%
  rownames()
column_order = match(column_order, rownames(subsobj@meta.data))

# Heatmap
ha_top = HeatmapAnnotation(sample_type = subsobj$sample_type,
                           sample_identifier = subsobj$sample_identifier,
                           col = list(sample_type = list_colors[["sample_type"]],
                                      sample_identifier = list_colors[["sample_identifier"]]))

# Heatmap
ht = Heatmap(as.matrix(mat_expr),
             heatmap_legend_param = list(title = "Expression", at = c(0, 1), 
                                         labels = c("low", "high")),
             col = list_colors[["expression"]],
             top_annotation = ha_top,
             show_column_names = FALSE,
             column_order = column_order,
             column_gap = unit(2, "mm"),
             cluster_rows = FALSE,
             row_title = NULL,
             row_names_gp = grid::gpar(fontsize = 14, fontface = "plain"),
             use_raster = FALSE,
             show_heatmap_legend = TRUE,
             border = TRUE)

ComplexHeatmap::draw(ht,
                     merge_legend = TRUE,
                     heatmap_legend_side = "bottom",
                     annotation_legend_side = "bottom")

HFSC

Settings

We load the HFSCs dataset :

sobj_hfsc = readRDS(paste0(data_dir, "/4_zoom/2_zoom_hfsc/hfsc_sobj.rds"))
sobj_hfsc

## An object of class Seurat 
## 15384 features across 1454 samples within 1 assay 
## Active assay: RNA (15384 features, 2000 variable features)
##  6 dimensional reductions calculated: RNA_pca, RNA_pca_24_tsne, RNA_pca_24_umap, harmony, harmony_24_umap, harmony_24_tsne

This is the projection name to visualize cells :

name2D = "harmony_24_tsne"

To represent results from differential expression, we load the analyses results :

list_results = readRDS(paste0(data_dir, "/4_zoom/2_zoom_hfsc/hfsc_list_results.rds"))

lapply(list_results, FUN = names)

## $cluster_0_8
## [1] "p_val"     "avg_logFC" "pct.1"     "pct.2"     "p_val_adj"
## 
## $cluster_2
## [1] "p_val"     "avg_logFC" "pct.1"     "pct.2"     "p_val_adj"
## 
## $cluster_1
## [1] "p_val"     "avg_logFC" "pct.1"     "pct.2"     "p_val_adj"
## 
## $cluster_3
## [1] "p_val"     "avg_logFC" "pct.1"     "pct.2"     "p_val_adj"

Figures

HFSCs on the full atlas :

sobj$is_hfsc = (colnames(sobj) %in% colnames(sobj_hfsc))

Seurat::DimPlot(sobj, reduction = name2D_atlas, pt.size = 0.000001,
                group.by = "is_hfsc", order = "TRUE") +
  ggplot2::scale_color_manual(values = c(color_markers[["HFSC"]], bg_color),
                              breaks = c(TRUE, FALSE)) +
  ggplot2::labs(title = "HFSCs",
                subtitle = paste0(ncol(sobj_hfsc), " cells")) +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_text(hjust = 0.5, face = "bold"),
                 plot.subtitle = element_text(hjust = 0.5)) +
  Seurat::NoAxes() + Seurat::NoLegend()

KRT15 expression :

Seurat::FeaturePlot(sobj, reduction = name2D_atlas, pt.size = 0.000001,
                    features = "KRT15") +
  ggplot2::scale_color_gradientn(colors = aquarius::color_gene) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() + Seurat::NoLegend()

Genes of interest :

genes = c("TGFB2", "ANGPTL7", "FGF18", "MGP", "EPCAM", "KRT75", "NOTCH3", "PTHLH")

plot_list = lapply(c(1:length(genes)), FUN = function(gene_id) {
  gene = genes[[gene_id]]
  
  sobj_hfsc$my_gene = Seurat::FetchData(sobj_hfsc, gene)[, 1] %>%
    aquarius::run_rescale(., new_min = 0, new_max = 10)
  
  Seurat::FeaturePlot(sobj_hfsc, features = "my_gene", reduction = name2D) +
    ggplot2::scale_color_gradientn(colors = aquarius::color_gene,
                                   breaks = seq(0, 10, by = 2.5),
                                   labels = c("min", rep("", 3), "max")) +
    ggplot2::labs(title = gene) +
    ggplot2::theme(aspect.ratio = 1,
                   plot.title = element_text(hjust = 0.5, size = 17),
                   plot.subtitle = element_text(hjust = 0.5, size = 15),
                   legend.text = element_text(size = 15),
                   legend.position = "none") +
    Seurat::NoAxes()
})

plot_list

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

plot_list = lapply(c("IFITM3", "MIF", "DDIT4"), FUN = function(one_gene) {
  p = aquarius::plot_split_dimred(sobj_hfsc,
                                  reduction = name2D,
                                  split_by = "sample_type",
                                  color_by = one_gene,
                                  color_palette = viridis::viridis(n = 10),
                                  main_pt_size = 0.6,
                                  bg_pt_size = 0.6,
                                  order = TRUE,
                                  bg_color = "gray95")
  p = patchwork::wrap_plots(p, nrow = 1) +
    patchwork::plot_layout(guides = "collect") +
    ggplot2::theme(legend.position = "right") &
    ggplot2::theme(plot.subtitle = element_blank())
  return(p)
})

plot_list

## [[1]]

## 
## [[2]]

## 
## [[3]]

Project name :

# Random order
set.seed(1234)
rnd_order = sample(colnames(sobj_hfsc), replace = FALSE, size = ncol(sobj_hfsc))

# Extract coordinates
cells_coord = sobj_hfsc@reductions[[name2D]]@cell.embeddings %>%
  as.data.frame() %>%
  `colnames<-`(c("Dim1", "Dim2"))
cells_coord$project_name = sobj_hfsc$project_name
cells_coord = cells_coord[(rnd_order), ]

# Plot
ggplot2::ggplot(cells_coord, aes(x = Dim1, y = Dim2, col = project_name)) +
  ggplot2::geom_point(size = 1.2) +
  ggplot2::scale_color_manual(values = sample_info$color,
                              breaks = sample_info$project_name) +
  ggplot2::theme_void() +
  ggplot2::theme(aspect.ratio = 1,
                 legend.position = "none")

Cluster :

Seurat::DimPlot(sobj_hfsc, reduction = name2D, pt.size = 1,
                group.by = "seurat_clusters", label = TRUE) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() +
  Seurat::NoLegend()

Heatmap with proportions :

cluster_markers = c("TGFB2", "ANGPTL7", "EPCAM", "KRT75", "NOTCH3", "PTHLH",
                    "percent.mt", "log_nCount_RNA")

## Bottom annotation : gene expression by cluster
ht_annot = Seurat::FetchData(sobj_hfsc, slot = "data", vars = cluster_markers) %>%
  as.data.frame()
ht_annot$clusters = sobj_hfsc$seurat_clusters
ht_annot = ht_annot %>%
  dplyr::group_by(clusters) %>%
  dplyr::summarise_all(funs('mean' = mean)) %>%
  as.data.frame() %>%
  dplyr::select(-clusters) %>%
  `colnames<-`(c(cluster_markers))

ha_bottom = ComplexHeatmap::HeatmapAnnotation(df = ht_annot,
                                              which = "column",
                                              show_legend = TRUE,
                                              col = setNames(nm = cluster_markers,
                                                             lapply(cluster_markers, FUN = color_fun)),
                                              annotation_name_side = "left")

## Right annotation : number of cells by dataset
ht_annot = table(sobj_hfsc$sample_identifier) %>%
  as.data.frame.table() %>%
  `colnames<-`(c("sample_identifier", "nb_cells")) %>%
  `rownames<-`(.$sample_identifier) %>%
  dplyr::select(-sample_identifier)

ha_right = ComplexHeatmap::HeatmapAnnotation(
  df = ht_annot,
  which = "row",
  show_legend = TRUE,
  annotation_name_side = "bottom",
  col = list(nb_cells  = circlize::colorRamp2(colors = RColorBrewer::brewer.pal(name = "Greys", n = 9),
                                              breaks = seq(from = range(ht_annot$nb_cells)[1],
                                                           to = range(ht_annot$nb_cells)[2],
                                                           length.out = 9))))

## Heatmap
ht = aquarius::plot_prop_heatmap(df = sobj_hfsc@meta.data[, c("sample_identifier", "seurat_clusters")],
                                 bottom_annotation = ha_bottom,
                                 right_annotation = ha_right,
                                 cluster_rows = TRUE,
                                 column_names_centered = TRUE,
                                 prop_margin = 1,
                                 row_names_gp = grid::gpar(names = sample_info$sample_identifier,
                                                           col = sample_info$color,
                                                           fontface = "bold"),
                                 row_title = "Sample",
                                 column_title = "Cluster")

ComplexHeatmap::draw(ht,
                     merge_legend = TRUE,
                     heatmap_legend_side = "bottom")

Dotplot of DE genes in clusters 0 and 8, between HS and HD :

# features_oi = rownames(list_results$cluster_0_8)
# features_oi = features_oi[!grepl(features_oi, pattern = "^RP")]

subsobj = subset(sobj_hfsc, seurat_clusters %in% c(0,8))

features_oi = c("HLA-C", "HLA-B", "IFITM3", "MIF", "JUNB", "MARCKSL1",
                "LTBP2", "PDK4", "DDIT4", "S100A11")

Seurat::DotPlot(subsobj,
                assay = "RNA",
                features = features_oi,
                group.by = "sample_type",
                scale = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::scale_color_gradientn(colors = aquarius::color_gene) +
  ggplot2::theme(axis.title.x = element_blank(),
                 axis.title.y = element_blank(),
                 axis.text.x = element_text(angle = 30, hjust = 1))

Heatmap for clulster 0 and 8 :

# Matrix
mat_expr = Seurat::GetAssayData(subsobj)
mat_expr = mat_expr[features_oi, ]
mat_expr = Matrix::t(mat_expr)
mat_expr = dynutils::scale_quantile(mat_expr) # between 0 and 1
mat_expr = Matrix::t(mat_expr)
dim(mat_expr) # genes x cells

## [1]  10 631

## Colors
list_colors = list()

# Heatmap
list_colors[["expression"]] = rev(RColorBrewer::brewer.pal(name = "RdBu", n = 9))

# Sample annotation (top annotation)
list_colors[["sample_type"]] = sample_type_colors
list_colors[["sample_identifier"]] = setNames(nm = sample_info$sample_identifier,
                                              sample_info$color)
# list_colors[["seurat_clusters"]] = setNames(aquarius::gg_color_hue(length(levels(subsobj$seurat_clusters))),
#                                             nm = levels(subsobj$seurat_clusters))
# Cells order
column_order = subsobj@meta.data %>%
  dplyr::arrange(sample_type, sample_identifier) %>%
  rownames()
column_order = match(column_order, rownames(subsobj@meta.data))

# Heatmap
ha_top = HeatmapAnnotation(sample_type = subsobj$sample_type,
                           sample_identifier = subsobj$sample_identifier,
                           # clusters = subsobj$seurat_clusters,
                           col = list(sample_type = list_colors[["sample_type"]],
                                      sample_identifier = list_colors[["sample_identifier"]]
                                      # clusters = list_colors[["seurat_clusters"]]
                           ))

# Heatmap
ht = Heatmap(as.matrix(mat_expr),
             heatmap_legend_param = list(title = "Expression", at = c(0, 1), 
                                         labels = c("low", "high")),
             col = list_colors[["expression"]],
             top_annotation = ha_top,
             show_column_names = FALSE,
             column_order = column_order,
             column_gap = unit(2, "mm"),
             cluster_rows = FALSE,
             row_title = NULL,
             row_names_gp = grid::gpar(fontsize = 14, fontface = "plain"),
             use_raster = FALSE,
             show_heatmap_legend = TRUE,
             border = TRUE)

ComplexHeatmap::draw(ht,
                     merge_legend = TRUE,
                     heatmap_legend_side = "bottom",
                     annotation_legend_side = "bottom")

Violin plot of IFITM3 :

table(subsobj$sample_type)

## 
##  HS  HD 
## 588  43

ifitm3_hs = subsobj@assays$RNA@data["IFITM3", subsobj$sample_type == "HS"]
ifitm3_hd = subsobj@assays$RNA@data["IFITM3", subsobj$sample_type == "HD"]
ifitm3_hs_VS_ifitm3_hd = stats::t.test(ifitm3_hs, ifitm3_hd)
ifitm3_hs_VS_ifitm3_hd

## 
##  Welch Two Sample t-test
## 
## data:  ifitm3_hs and ifitm3_hd
## t = 7.0204, df = 47.675, p-value = 7.081e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.5509680 0.9933324
## sample estimates:
## mean of x mean of y 
##  1.775647  1.003497

Seurat::VlnPlot(subsobj, group.by = "sample_type",
                features = "IFITM3", cols = sample_type_colors) +
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

Violin plot of PDK4 :

table(subsobj$sample_type)

## 
##  HS  HD 
## 588  43

pdk4_hs = subsobj@assays$RNA@data["PDK4", subsobj$sample_type == "HS"]
pdk4_hd = subsobj@assays$RNA@data["PDK4", subsobj$sample_type == "HD"]
pdk4_hs_VS_pdk4_hd = stats::t.test(pdk4_hs, pdk4_hd)
pdk4_hs_VS_pdk4_hd

## 
##  Welch Two Sample t-test
## 
## data:  pdk4_hs and pdk4_hd
## t = -10.934, df = 43.02, p-value = 5.329e-14
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.523546 -1.737614
## sample estimates:
## mean of x mean of y 
## 0.1615592 2.2921393

Seurat::VlnPlot(subsobj, group.by = "sample_type",
                features = "PDK4", cols = sample_type_colors) +
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

IBL and ORS

Settings

We load the IBL + ORS dataset :

sobj_iblors = readRDS(paste0(data_dir, "/4_zoom/3_zoom_iblmors/iblmors_sobj.rds"))
sobj_iblors

## An object of class Seurat 
## 16701 features across 3532 samples within 1 assay 
## Active assay: RNA (16701 features, 2000 variable features)
##  6 dimensional reductions calculated: RNA_pca, RNA_pca_20_tsne, RNA_pca_20_umap, harmony, harmony_20_umap, harmony_20_tsne

This is the projection name to visualize cells :

name2D = "harmony_20_tsne"

To represent results from differential expression, we load the analyses results :

list_results = readRDS(paste0(data_dir, "/4_zoom/3_zoom_iblmors/iblmors_list_results.rds"))

lapply(list_results, FUN = names)

## $IBL_vs_ORS
## [1] "mark"        "enrichr_ibl" "enrichr_ors" "gsea"       
## 
## $cluster5_vs_ORS
## [1] "mark"       "enrichr_up" "enrichr_dn" "gsea"      
## 
## $IBL_HS_vs_HD
## [1] "mark"       "enrichr_hs" "enrichr_hd" "gsea"      
## 
## $ORS_HS_vs_HD
## [1] "mark"       "enrichr_hs" "enrichr_hd" "gsea"

Preparation

We defined cluster type :

cluster_type = table(sobj_iblors$cell_type, sobj_iblors$seurat_clusters) %>%
  prop.table(., margin = 2) %>%
  apply(., 2, which.max)
cluster_type = setNames(nm = names(cluster_type),
                        levels(sobj_iblors$cell_type)[cluster_type])

sobj_iblors$cluster_type = cluster_type[sobj_iblors$seurat_clusters]
sobj_iblors$cluster_type = factor(sobj_iblors$cluster_type,
                                  levels = levels(sobj_iblors$cell_type))

Figures

IBL + ORS on the full atlas :

sobj$cell_bc = colnames(sobj)
sobj_iblors$cell_bc = colnames(sobj_iblors)
sobj$is_iblors = dplyr::left_join(sobj@meta.data[, c("cell_bc", "percent.mt")],
                                  sobj_iblors@meta.data[, c("cell_bc", "cluster_type")],
                                  by = "cell_bc")[, "cluster_type"]

Seurat::DimPlot(sobj, reduction = name2D_atlas, pt.size = 0.000001,
                group.by = "is_iblors", order = "TRUE") +
  ggplot2::scale_color_manual(values = c(color_markers[c("IBL", "ORS")], bg_color),
                              breaks = c("IBL", "ORS", NA), na.value = bg_color) +
  ggplot2::labs(title = "IBL + ORS",
                subtitle = paste0(ncol(sobj_iblors), " cells")) +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_text(hjust = 0.5, face = "bold"),
                 plot.subtitle = element_text(hjust = 0.5)) +
  Seurat::NoAxes() + Seurat::NoLegend()

Project name :

# Random order
set.seed(1234)
rnd_order = sample(colnames(sobj_iblors), replace = FALSE, size = ncol(sobj_iblors))

# Extract coordinates
cells_coord = sobj_iblors@reductions[[name2D]]@cell.embeddings %>%
  as.data.frame() %>%
  `colnames<-`(c("Dim1", "Dim2"))
cells_coord$project_name = sobj_iblors$project_name
cells_coord = cells_coord[(rnd_order), ]

# Plot
ggplot2::ggplot(cells_coord, aes(x = Dim1, y = Dim2, col = project_name)) +
  ggplot2::geom_point(size = 1.2) +
  ggplot2::scale_color_manual(values = sample_info$color,
                              breaks = sample_info$project_name) +
  ggplot2::theme_void() +
  ggplot2::theme(aspect.ratio = 1,
                 legend.position = "none")

Cluster :

Seurat::DimPlot(sobj_iblors, reduction = name2D, pt.size = 1,
                group.by = "seurat_clusters", label = TRUE) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() +
  Seurat::NoLegend()

Cluster type :

Seurat::DimPlot(sobj_iblors, reduction = name2D, pt.size = 1,
                group.by = "cluster_type", cols = color_markers) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() +
  Seurat::NoLegend()

Cluster type split by sample type :

plot_list = aquarius::plot_split_dimred(sobj_iblors, reduction = name2D,
                                        group_by = "cluster_type",
                                        group_color = color_markers,
                                        split_by = "sample_type",
                                        bg_pt_size = 1, main_pt_size = 1,
                                        bg_color = bg_color)

patchwork::wrap_plots(plot_list) &
  Seurat::NoLegend()

DE genes between IBL and ORS :

mark = list_results$IBL_vs_ORS$mark
mark$gene_name = rownames(mark)
mark_label = rbind(
  # up-regulated in IBL
  mark %>% dplyr::top_n(., n = 20, wt = avg_logFC),
  # up-regulated in ORS
  mark %>% dplyr::top_n(., n = 20, wt = -avg_logFC),
  # representative and selective for IBL
  mark %>% dplyr::top_n(., n = 20, wt = (pct.1 - pct.2)),
  # representative and selective for ORS
  mark %>% dplyr::top_n(., n = 20, wt = -(pct.1 - pct.2))) %>%
  dplyr::distinct()
mark_label = mark_label[!grepl(rownames(mark_label), pattern = "^MT"), ]

ggplot2::ggplot(mark, aes(x = pct.1, y = pct.2, col = avg_logFC)) +
  ggplot2::geom_abline(slope = 1, intercept = 0, lty = 2) +
  ggplot2::geom_point() +
  ggrepel::geom_label_repel(data = mark_label, max.overlaps = Inf,
                            aes(x = pct.1, y = pct.2, label = gene_name),
                            col = "black", fill = NA, size = 3, label.size = NA) +
  ggplot2::labs(title = "Differentially expressed genes",
                subtitle = "between IBL and ORS") +
  ggplot2::scale_color_gradient2(low = aquarius::color_cnv[1],
                                 mid = aquarius::color_cnv[2],
                                 high = aquarius::color_cnv[3],
                                 midpoint = 0) +
  ggplot2::theme_classic() +
  ggplot2::theme(plot.title = element_text(hjust = 0.5),
                 plot.subtitle = element_text(hjust = 0.5))

GSEA for all genes, in IBL, between HS and HD :

# gs_oi = c("REACTOME_EXPORT_OF_VIRAL_RIBONUCLEOPROTEINS_FROM_NUCLEUS",
#           "GOBP_RESPONSE_TO_UV_C",
#           "GOBP_POSITIVE_REGULATION_OF_HISTONE_H3_K9_METHYLATION",
#           "GOBP_HISTONE_H3_K14_ACETYLATION",
#           "GOBP_MITOPHAGY")

list_results$IBL_HS_vs_HD$gsea@result %>%
  # dplyr::filter(ID %in% gs_oi) %>%
  dplyr::filter(pvalue <  0.05) %>%
  dplyr::top_n(., n = 200, wt = abs(NES)) %>%
  dplyr::mutate(too_long = ifelse(nchar(ID) > 60, yes = TRUE, no = FALSE)) %>%
  dplyr::mutate(ID = stringr::str_sub(ID, end = 60)) %>%
  dplyr::mutate(ID = ifelse(too_long, yes = paste0(ID, "..."), no = ID)) %>%
  aquarius::gsea_plot(show_legend = TRUE) +
  ggplot2::labs(title = "GSEA using all genes (count matrix)") +
  ggplot2::theme(plot.title = element_text(size = 20))

GSEA for all genes, in ORS, between HS and HD :

list_results$ORS_HS_vs_HD$gsea@result %>%
  # dplyr::filter(ID %in% gs_oi) %>%
  dplyr::filter(pvalue <  0.05) %>%
  dplyr::top_n(., n = 200, wt = abs(NES)) %>%
  dplyr::mutate(too_long = ifelse(nchar(ID) > 70, yes = TRUE, no = FALSE)) %>%
  dplyr::mutate(ID = stringr::str_sub(ID, end = 70)) %>%
  dplyr::mutate(ID = ifelse(too_long, yes = paste0(ID, "..."), no = ID)) %>%
  aquarius::gsea_plot(show_legend = TRUE) +
  ggplot2::labs(title = "GSEA using all genes (count matrix)") +
  ggplot2::theme(plot.title = element_text(size = 20))

Violin plot for IBL :

subsobj = subset(sobj_iblors, cluster_type == "IBL")

Seurat::VlnPlot(subsobj, group.by = "sample_type",
                features = c("DUSP1", "DDIT4", "MIF", "LGALS7", "ARF5", "S100A9"),
                cols = sample_type_colors, ncol = 3) &
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

Violin plot for ORS :

subsobj = subset(sobj_iblors, cluster_type == "ORS")

Seurat::VlnPlot(subsobj, group.by = "sample_type",
                features = c("DUSP1", "CLDN1", "KLF6", "CTGF", "CCL2", "S100A9", "IFI27", "IFITM3"),
                cols = sample_type_colors, ncol = 4) &
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

Split by sample :

Seurat::VlnPlot(subsobj, group.by = "sample_identifier",
                features = "IFI27", cols = sample_info$color) +
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = "none")

Dotplot of DE genes in ORS, between cluster 5 and other ORS :

subsobj = subset(sobj_iblors, cluster_type == "ORS")
subsobj = subset(subsobj, seurat_clusters != 7)

# features_oi = list_results$cluster5_vs_ORS$mark %>%
#   dplyr::filter(p_val_adj < 0.05) %>%
#   dplyr::filter(abs(avg_logFC) > 0.5) %>%
#   dplyr::arrange(avg_logFC) %>%
#   rownames()

# features_oi = c("YBX3", "TXNIP", "KRT15", "NEAT1", "COL17A1", "FXYD3", "KRT14",
#                 "CXCL14", "MT2A", "MT1E", "FOS", "JUNB", "AQP3", "GLUL", "KLF5",
#                 "GSTM3", "ALDH3A1", "HLA-C", "LGALS7B", "SLC38A2", "EHF", "KLF3",
#                 "IL18", "CLEC2B", "IL20RB", "GPSM2", "DAAM1", "DUSP1", "KLF6",
#                 "ZFP36", "ATF3", "RHOB", "MT1X", "KLF4", "ETS2", "NFKBIZ", "ID1",
#                 "RNASET2", "HOPX", "WNT4", "POU3F1", "SPRY1", "AR", "THSD4", "PDGFC",
#                 "KRT31", "WFDC2", "SPINK5", "SLPI", "IL1R2", "PLAT", "TSC22D3",
#                 "KLF9", "WNT3", "FGFR3", "LAMB4", "IFI27", "DCN", "LY6D", "IGFBP3", "WFDC5",
#                 
#                 "S100A14", "TMSB10", "ACTG1", "ACTB", "CSTB", "APOE", "CALD1", "SOX4",
#                 "STMN1", "LMO4", "CEBPB", "TMEM45A", "CTSB", "GPX2", "C1QTNF12", "GJB6",
#                 "RBP1", "TPM1", "KRT17", "CALML3", "PTN", "DAPK2",
#                 "EGLN3", "FILIP1L", "FOXC1", "ADGRL3", "FST", "EFNB2", "SEMA5A",
#                 "FGFR1", "DACH1", "EGR2", "CLDN1", "CTGF", "DEFB1", "CARD18", "KRT6A",
#                 "S100A7", "KRT6B", "MGST1", "CHI3L1", "CCL19", "LGALS1", "CYP1B1", "VIM")

features_oi = c("YBX3", "TXNIP", "KRT14", "KRT15", "NEAT1",
                "FXYD3", "MT2A", "MT1E", "MT1X", "AQP3", "GLUL",
                # "HALLMARK_TNFA_SIGNALING_VIA_NFKB"
                "FOS", "JUNB", "DUSP1", "ZFP36", "NFKBIZ",
                "ATF3", "RHOB",  "ETS2", "IL18", "KLF4", "KLF6", "KLF9",
                "KLF3", "KLF5", "COL17A1", "THSD4", "WNT3", "WNT4", "SLPI", "PLAT",
                "LAMB4", "DCN", "SPINK5",
                "GSTM3", "ALDH3A1",  "LGALS7B", "SLC38A2", "EHF",  "CLEC2B",
                "IL20RB", "IL1R2", "IFI27", "CXCL14", "HLA-C", "GPSM2", "DAAM1",   "ID1",
                "RNASET2", "HOPX", "POU3F1", "SPRY1", "AR", "PDGFC",
                "WFDC2", "WFDC5", "TSC22D3", "FGFR3",  "LY6D", "IGFBP3", 
                
                # Other ORS
                "APOE", "CTSB", "CALD1", "SOX4",
                "STMN1", "LMO4", "CEBPB", "TMEM45A", "GPX2", "C1QTNF12", "GJB6",
                "KRT6A", "KRT17", "RBP1", "CALML3", "PTN", "DAPK2",
                "EGLN3", "FILIP1L", "ADGRL3", "FST", "EFNB2", "SEMA5A",
                "FGFR1", "EGR2", "CLDN1", "DEFB1", "CARD18", "MGST1")

Seurat::DotPlot(subsobj,
                assay = "RNA",
                features = features_oi,
                group.by = "sample_type",
                scale = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::scale_color_gradientn(colors = aquarius::color_gene) +
  ggplot2::theme(axis.title.x = element_blank(),
                 axis.title.y = element_blank(),
                 axis.text.x = element_text(angle = 30, hjust = 1))

Heatmap for cluster 5 vs other ORS :

# Matrix
mat_expr = Seurat::GetAssayData(subsobj)
mat_expr = mat_expr[features_oi, ]
mat_expr = Matrix::t(mat_expr)
mat_expr = dynutils::scale_quantile(mat_expr) # between 0 and 1
mat_expr = Matrix::t(mat_expr)
dim(mat_expr) # genes x cells

## [1]   89 1793

## Colors
list_colors = list()
list_colors[["expression"]] = rev(RColorBrewer::brewer.pal(name = "RdBu", n = 9))
list_colors[["sample_type"]] = sample_type_colors
list_colors[["sample_identifier"]] = setNames(nm = sample_info$sample_identifier,
                                              sample_info$color)
list_colors[["seurat_clusters"]] = setNames(aquarius::gg_color_hue(length(levels(subsobj$seurat_clusters))),
                                            nm = levels(subsobj$seurat_clusters))
list_colors[["seurat_clusters"]] = list_colors[["seurat_clusters"]][unique(subsobj$seurat_clusters)]

# Cells order
column_order = subsobj@meta.data %>%
  dplyr::arrange(sample_type, sample_identifier) %>%
  rownames()
column_order = match(column_order, rownames(subsobj@meta.data))

# Annotation
ha_top = HeatmapAnnotation(sample_type = subsobj$sample_type,
                           sample_identifier = subsobj$sample_identifier,
                           clusters = subsobj$seurat_clusters,
                           col = list(sample_type = list_colors[["sample_type"]],
                                      sample_identifier = list_colors[["sample_identifier"]],
                                      clusters = list_colors[["seurat_clusters"]]))

# Heatmap
ht = Heatmap(as.matrix(mat_expr),
             heatmap_legend_param = list(title = "Expression", at = c(0, 1), 
                                         labels = c("low", "high")),
             col = list_colors[["expression"]],
             top_annotation = ha_top,
             # Cell grouping
             column_split = subsobj$sample_type %>% as.character(),
             cluster_columns = TRUE,
             column_title = NULL,
             show_column_dend = FALSE,
             show_column_names = FALSE,
             # Genes
             cluster_rows = FALSE,
             row_names_gp = grid::gpar(fontsize = 14, fontface = "plain"),
             # Style
             use_raster = FALSE,
             show_heatmap_legend = TRUE,
             border = TRUE)

ComplexHeatmap::draw(ht,
                     merge_legend = TRUE,
                     heatmap_legend_side = "bottom",
                     annotation_legend_side = "bottom")

Genes of interest :

genes = c("KRT16", "COL17A1", "DST", "KRT6B", "IL1R2", "WNT3",
          "IFI27", "CXCL14", "IGFBP3", "KRT15", "CD200")

plot_list = lapply(c(1:length(genes)), FUN = function(gene_id) {
  gene = genes[[gene_id]]
  
  sobj_iblors$my_gene = Seurat::FetchData(sobj_iblors, gene)[, 1] %>%
    aquarius::run_rescale(., new_min = 0, new_max = 10)
  
  Seurat::FeaturePlot(sobj_iblors, features = "my_gene", reduction = name2D, pt.size = 0.25) +
    ggplot2::scale_color_gradientn(colors = aquarius::color_gene,
                                   breaks = seq(0, 10, by = 2.5),
                                   labels = c("min", rep("", 3), "max")) +
    ggplot2::labs(title = gene) +
    ggplot2::theme(aspect.ratio = 1,
                   plot.title = element_text(hjust = 0.5, size = 17),
                   plot.subtitle = element_text(hjust = 0.5, size = 15),
                   legend.text = element_text(size = 15),
                   legend.position = "none") +
    Seurat::NoAxes()
})

plot_list

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

HFSCs to IBL and ORS

Settings

We load the merged dataset :

sobj_traj = readRDS(paste0(data_dir, "/4_zoom/4_zoom_hfsc_iblmors/hfsc_iblmors_sobj_traj_tinga.rds"))
sobj_traj

## An object of class Seurat 
## 17050 features across 4986 samples within 1 assay 
## Active assay: RNA (17050 features, 2000 variable features)
##  8 dimensional reductions calculated: RNA_pca, RNA_pca_18_tsne, RNA_pca_18_umap, harmony, harmony_18_umap, harmony_18_tsne, harmony_dm, harmony_dm_5_umap

This is the projection name to visualize cells :

name2D = "harmony_dm"

We load the trajectory object for visualisation purpose :

my_traj = readRDS(paste0(data_dir, "/4_zoom/4_zoom_hfsc_iblmors/hfsc_iblmors_my_traj_tinga.rds"))
class(my_traj)

## [1] "dynwrap::with_dimred"     "dynwrap::with_trajectory"
## [3] "dynwrap::data_wrapper"    "list"

Preparation

We defined cell type based on individual object :

sobj_iblors$cell_bc = colnames(sobj_iblors)
sobj_traj$cell_bc = colnames(sobj_traj)
sobj_traj$cluster_type = dplyr::left_join(sobj_traj@meta.data[, c("cell_bc", "percent.mt")],
                                          sobj_iblors@meta.data[, c("cell_bc", "cluster_type")],
                                          by = "cell_bc")[, "cluster_type"] %>% as.character()
sobj_traj$cluster_type = ifelse(colnames(sobj_traj) %in% colnames(sobj_hfsc),
                                yes = "HFSC",
                                no = sobj_traj$cluster_type) %>%
  as.factor()

Figures

Cells on the full atlas :

sobj$cell_bc = colnames(sobj)
sobj_traj$cell_bc = colnames(sobj_traj)
sobj$is_traj = dplyr::left_join(sobj@meta.data[, c("cell_bc", "percent.mt")],
                                sobj_traj@meta.data[, c("cell_bc", "cluster_type")],
                                by = "cell_bc")[, "cluster_type"]

Seurat::DimPlot(sobj, reduction = name2D_atlas, pt.size = 0.000001,
                group.by = "is_traj", order = levels(sobj_traj$cluster_type)) +
  ggplot2::scale_color_manual(values = c(color_markers[c("IBL", "ORS", "HFSC")], bg_color),
                              breaks = c("IBL", "ORS", "HFSC", NA), na.value = bg_color) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() + Seurat::NoLegend()

Project name :

# Random order
set.seed(1234)
rnd_order = sample(colnames(sobj_traj), replace = FALSE, size = ncol(sobj_traj))

# Extract coordinates
cells_coord = sobj_traj@reductions[[name2D]]@cell.embeddings %>%
  as.data.frame() %>%
  `colnames<-`(c("Dim1", "Dim2"))
cells_coord$project_name = sobj_traj$project_name
cells_coord = cells_coord[(rnd_order), ]

# Plot
ggplot2::ggplot(cells_coord, aes(x = Dim1, y = Dim2, col = project_name)) +
  ggplot2::geom_point(size = 1.2) +
  ggplot2::scale_color_manual(values = sample_info$color,
                              breaks = sample_info$project_name) +
  ggplot2::theme_void() +
  ggplot2::theme(aspect.ratio = 1,
                 legend.position = "none")

Cluster :

Seurat::DimPlot(sobj_traj, reduction = name2D, pt.size = 1,
                group.by = "seurat_clusters", label = TRUE) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() +
  Seurat::NoLegend()

Cluster type :

Seurat::DimPlot(sobj_traj, reduction = name2D, pt.size = 1,
                group.by = "cluster_type", cols = color_markers) +
  ggplot2::theme(aspect.ratio = 1) +
  Seurat::NoAxes() +
  Seurat::NoLegend()

Pseudotime :

Seurat::FeaturePlot(sobj_traj, reduction = name2D, pt.size = 0.5,
                    features = "pseudotime") +
  ggplot2::scale_color_gradientn(colors = viridis::viridis(n = 100)) +
  ggplot2::lims(x = range(sobj_traj@reductions[[name2D]]@cell.embeddings[, 1]),
                y = range(sobj_traj@reductions[[name2D]]@cell.embeddings[, 2])) +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_blank()) +
  Seurat::NoAxes()

Pseudotime with dynplot’s function :

dynplot::plot_dimred(trajectory = my_traj,
                     dimred = sobj_traj[[name2D]]@cell.embeddings,
                     # Cells
                     color_cells = 'pseudotime',
                     size_cells = 1.6,
                     border_radius_percentage = 0,
                     # Trajectory
                     plot_trajectory = TRUE,
                     color_trajectory = "none",
                     label_milestones = FALSE,
                     size_milestones = 0,
                     size_transitions = 1)

R Session

show

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
## 
## locale:
## [1] C
## 
## attached base packages:
##  [1] parallel  stats4    grid      stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] org.Mm.eg.db_3.10.0   AnnotationDbi_1.48.0  IRanges_2.20.2       
##  [4] S4Vectors_0.24.4      Biobase_2.46.0        BiocGenerics_0.32.0  
##  [7] ComplexHeatmap_2.14.0 ggplot2_3.3.5         patchwork_1.1.2      
## [10] dplyr_1.0.7          
## 
## loaded via a namespace (and not attached):
##   [1] softImpute_1.4              graphlayouts_0.7.0         
##   [3] pbapply_1.4-2               lattice_0.20-41            
##   [5] haven_2.3.1                 dyndimred_1.0.3            
##   [7] vctrs_0.3.8                 usethis_2.0.1              
##   [9] dynwrap_1.2.1               blob_1.2.1                 
##  [11] survival_3.2-13             prodlim_2019.11.13         
##  [13] dynutils_1.0.5              later_1.3.0                
##  [15] DBI_1.1.1                   R.utils_2.11.0             
##  [17] SingleCellExperiment_1.8.0  rappdirs_0.3.3             
##  [19] uwot_0.1.8                  dqrng_0.2.1                
##  [21] jpeg_0.1-8.1                zlibbioc_1.32.0            
##  [23] pspline_1.0-18              pcaMethods_1.78.0          
##  [25] mvtnorm_1.1-1               htmlwidgets_1.5.4          
##  [27] GlobalOptions_0.1.2         future_1.22.1              
##  [29] UpSetR_1.4.0                laeken_0.5.2               
##  [31] leiden_0.3.3                clustree_0.4.3             
##  [33] lmds_0.1.0                  scater_1.14.6              
##  [35] irlba_2.3.3                 DEoptimR_1.0-9             
##  [37] tidygraph_1.1.2             Rcpp_1.0.9                 
##  [39] readr_2.0.2                 KernSmooth_2.23-17         
##  [41] carrier_0.1.0               promises_1.1.0             
##  [43] gdata_2.18.0                DelayedArray_0.12.3        
##  [45] limma_3.42.2                graph_1.64.0               
##  [47] RcppParallel_5.1.4          Hmisc_4.4-0                
##  [49] fs_1.5.2                    RSpectra_0.16-0            
##  [51] fastmatch_1.1-0             ranger_0.12.1              
##  [53] digest_0.6.25               png_0.1-7                  
##  [55] sctransform_0.2.1           cowplot_1.0.0              
##  [57] DOSE_3.12.0                 here_1.0.1                 
##  [59] TInGa_0.0.0.9000            dynplot_1.1.0              
##  [61] ggraph_2.0.3                pkgconfig_2.0.3            
##  [63] GO.db_3.10.0                DelayedMatrixStats_1.8.0   
##  [65] gower_0.2.1                 ggbeeswarm_0.6.0           
##  [67] iterators_1.0.12            DropletUtils_1.6.1         
##  [69] reticulate_1.26             clusterProfiler_3.14.3     
##  [71] SummarizedExperiment_1.16.1 circlize_0.4.15            
##  [73] beeswarm_0.4.0              GetoptLong_1.0.5           
##  [75] xfun_0.35                   bslib_0.3.1                
##  [77] zoo_1.8-10                  tidyselect_1.1.0           
##  [79] GA_3.2                      reshape2_1.4.4             
##  [81] purrr_0.3.4                 ica_1.0-2                  
##  [83] pcaPP_1.9-73                viridisLite_0.3.0          
##  [85] rtracklayer_1.46.0          rlang_1.0.2                
##  [87] hexbin_1.28.1               jquerylib_0.1.4            
##  [89] dyneval_0.9.9               glue_1.4.2                 
##  [91] waldo_0.3.1                 RColorBrewer_1.1-2         
##  [93] matrixStats_0.56.0          stringr_1.4.0              
##  [95] lava_1.6.7                  europepmc_0.3              
##  [97] DESeq2_1.26.0               recipes_0.1.17             
##  [99] labeling_0.3                httpuv_1.5.2               
## [101] class_7.3-17                BiocNeighbors_1.4.2        
## [103] DO.db_2.9                   annotate_1.64.0            
## [105] jsonlite_1.7.2              XVector_0.26.0             
## [107] bit_4.0.4                   mime_0.9                   
## [109] aquarius_0.1.5              Rsamtools_2.2.3            
## [111] gridExtra_2.3               gplots_3.0.3               
## [113] stringi_1.4.6               processx_3.5.2             
## [115] gsl_2.1-6                   bitops_1.0-6               
## [117] cli_3.0.1                   batchelor_1.2.4            
## [119] RSQLite_2.2.0               randomForest_4.6-14        
## [121] tidyr_1.1.4                 data.table_1.14.2          
## [123] rstudioapi_0.13             GenomicAlignments_1.22.1   
## [125] nlme_3.1-147                qvalue_2.18.0              
## [127] scran_1.14.6                locfit_1.5-9.4             
## [129] scDblFinder_1.1.8           listenv_0.8.0              
## [131] ggthemes_4.2.4              gridGraphics_0.5-0         
## [133] R.oo_1.24.0                 dbplyr_1.4.4               
## [135] TTR_0.24.2                  readxl_1.3.1               
## [137] lifecycle_1.0.1             timeDate_3043.102          
## [139] ggpattern_0.3.1             munsell_0.5.0              
## [141] cellranger_1.1.0            R.methodsS3_1.8.1          
## [143] proxyC_0.1.5                visNetwork_2.0.9           
## [145] caTools_1.18.0              codetools_0.2-16           
## [147] GenomeInfoDb_1.22.1         vipor_0.4.5                
## [149] lmtest_0.9-38               msigdbr_7.5.1              
## [151] htmlTable_1.13.3            triebeard_0.3.0            
## [153] lsei_1.2-0                  xtable_1.8-4               
## [155] ROCR_1.0-7                  BiocManager_1.30.10        
## [157] scatterplot3d_0.3-41        abind_1.4-5                
## [159] farver_2.0.3                parallelly_1.28.1          
## [161] RANN_2.6.1                  askpass_1.1                
## [163] GenomicRanges_1.38.0        RcppAnnoy_0.0.16           
## [165] tibble_3.1.5                ggdendro_0.1-20            
## [167] cluster_2.1.0               future.apply_1.5.0         
## [169] Seurat_3.1.5                dendextend_1.15.1          
## [171] Matrix_1.3-2                ellipsis_0.3.2             
## [173] prettyunits_1.1.1           lubridate_1.7.9            
## [175] ggridges_0.5.2              igraph_1.2.5               
## [177] RcppEigen_0.3.3.7.0         fgsea_1.12.0               
## [179] remotes_2.4.2               scBFA_1.0.0                
## [181] destiny_3.0.1               VIM_6.1.1                  
## [183] testthat_3.1.0              htmltools_0.5.2            
## [185] BiocFileCache_1.10.2        yaml_2.2.1                 
## [187] utf8_1.1.4                  plotly_4.9.2.1             
## [189] XML_3.99-0.3                ModelMetrics_1.2.2.2       
## [191] e1071_1.7-3                 foreign_0.8-76             
## [193] withr_2.5.0                 fitdistrplus_1.0-14        
## [195] BiocParallel_1.20.1         xgboost_1.4.1.1            
## [197] bit64_4.0.5                 foreach_1.5.0              
## [199] robustbase_0.93-9           Biostrings_2.54.0          
## [201] GOSemSim_2.13.1             rsvd_1.0.3                 
## [203] memoise_2.0.0               evaluate_0.18              
## [205] forcats_0.5.0               rio_0.5.16                 
## [207] geneplotter_1.64.0          tzdb_0.1.2                 
## [209] caret_6.0-86                ps_1.6.0                   
## [211] DiagrammeR_1.0.6.1          curl_4.3                   
## [213] fdrtool_1.2.15              fansi_0.4.1                
## [215] highr_0.8                   urltools_1.7.3             
## [217] xts_0.12.1                  GSEABase_1.48.0            
## [219] acepack_1.4.1               edgeR_3.28.1               
## [221] checkmate_2.0.0             scds_1.2.0                 
## [223] cachem_1.0.6                npsurv_0.4-0               
## [225] babelgene_22.3              rjson_0.2.20               
## [227] openxlsx_4.1.5              ggrepel_0.9.1              
## [229] clue_0.3-60                 rprojroot_2.0.2            
## [231] stabledist_0.7-1            tools_3.6.3                
## [233] sass_0.4.0                  nichenetr_1.1.1            
## [235] magrittr_2.0.1              RCurl_1.98-1.2             
## [237] proxy_0.4-24                car_3.0-11                 
## [239] ape_5.3                     ggplotify_0.0.5            
## [241] xml2_1.3.2                  httr_1.4.2                 
## [243] assertthat_0.2.1            rmarkdown_2.18             
## [245] boot_1.3-25                 globals_0.14.0             
## [247] R6_2.4.1                    Rhdf5lib_1.8.0             
## [249] nnet_7.3-14                 RcppHNSW_0.2.0             
## [251] progress_1.2.2              genefilter_1.68.0          
## [253] statmod_1.4.34              gtools_3.8.2               
## [255] shape_1.4.6                 HDF5Array_1.14.4           
## [257] BiocSingular_1.2.2          rhdf5_2.30.1               
## [259] splines_3.6.3               AUCell_1.8.0               
## [261] carData_3.0-4               colorspace_1.4-1           
## [263] generics_0.1.0              base64enc_0.1-3            
## [265] dynfeature_1.0.0            smoother_1.1               
## [267] gridtext_0.1.1              pillar_1.6.3               
## [269] tweenr_1.0.1                sp_1.4-1                   
## [271] ggplot.multistats_1.0.0     rvcheck_0.1.8              
## [273] GenomeInfoDbData_1.2.2      plyr_1.8.6                 
## [275] gtable_0.3.0                zip_2.2.0                  
## [277] knitr_1.41                  latticeExtra_0.6-29        
## [279] biomaRt_2.42.1              fastmap_1.1.0              
## [281] ADGofTest_0.3               copula_1.0-0               
## [283] doParallel_1.0.15           vcd_1.4-8                  
## [285] babelwhale_1.0.1            openssl_1.4.1              
## [287] scales_1.1.1                backports_1.2.1            
## [289] ipred_0.9-12                enrichplot_1.6.1           
## [291] hms_1.1.1                   ggforce_0.3.1              
## [293] Rtsne_0.15                  shiny_1.7.1                
## [295] numDeriv_2016.8-1.1         polyclip_1.10-0            
## [297] lazyeval_0.2.2              Formula_1.2-3              
## [299] tsne_0.1-3                  crayon_1.3.4               
## [301] MASS_7.3-54                 pROC_1.16.2                
## [303] viridis_0.5.1               dynparam_1.0.0             
## [305] rpart_4.1-15                zinbwave_1.8.0             
## [307] compiler_3.6.3              ggtext_0.1.0